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It is generally believed that the inhomogeneous Larkin-Ovchinnikov-Fulde-Ferrell (LOFF) phase appears in 
a color superconductor when the pairing between different quark flavors is under the circumstances of mis¬ 
matched Fermi surfaces. However, the real crystal structure of the LOFF phase is still unclear because an exact 
treatment of 3D crystal structures is rather difficult. In this work we present a solid-state-like calculation of 
the ground-state energy of the hody-centered cubic (BCC) structure for two-flavor pairing by diagonalizing the 
Hamiltonian matrix in the Bloch space without assuming a small amplitude of the order parameter. We develop 
a computational scheme to overcome the difficulties in diagonalizing huge matrices. Our results show that the 
BCC structure is energetically more favorable than the ID modulation in a narrow window around the con¬ 
ventional LOFF-normal phase transition point, which indicates the significance of the higher-order terms in the 
Ginzburg-Landau approach. 

PACS numbers: 12.38.-t, 21.65.Qr, 74.20.Fg, 03.75.Hh 


I. INTRODUCTION 

The ground state of exotic fermion Cooper pairing with 
mismatched Fermi surfaces is a longstanding problem in the 
theory of superconductivity iQl]. In electronic superconduc¬ 
tors, the mismatched Fermi surfaces are normally induced by 
the Zeeman energy splitting 25p in a magnetic field. For s- 
wave pairing at weak coupling, it is known that, at a critical 
field dpi - 0.707Ao where Aq is the pairing gap at vanish¬ 
ing mismatch, a first-order phase transition from the gapped 
BCS state to the normal state occurs i]. Further theoretical 
studies showed that the inhomogeneous Larkin-Ovchinnikov- 
Fulde-Ferrell (LOFF) state can survive in a narrow window 
< 6 p < 6 p 2 , where the upper critical field 6 p 2 - 0.754Ao 
13, 3]. However, since the thermodynamic critical field is 
much lower than dpi due to strong orbit effect, it is rather hard 
to observe the LOFF state in ordinary superconductors iQ]. In 
recent years, experimental evidences for the LOFF state in 
some superconducting materials have been reported M. 

On the other hand, exotic pairing phases have promoted 
new interest in the studies of densequark matter under the 
circumstances of compact stars H-tlMl and ultracold atomic 
Fermi gases with population imbalance Color super¬ 

conductivity in dense quark matter appears due to the attrac¬ 
tive interactions in certain diquark channels lE^ - E^ . Because 
of the constraints from Beta equilibrium and electric charge 
neutrality, different quark flavors (u, d, and i) acquire mis¬ 
matched Fermi surfaces. Quark color superconductors under 
compact-star constraints as well as atomic Fermi gases with 
population imbalance therefore provide rather clean systems 
to realize the long-sought exotic LOFF phase. 

Around the tricritical point in the temperature-mismatch 
phase diagram, the LOFF phase can be studied rigorously by 
using the Ginzburg-Laudau (GL) analysis since both the rap 
parameter and the pair momentum are vanishingly small 111]. 
It was found that the solution with two antipodal wave vectors 
is the preferred one . However, the real ground state of 

the LOFF phase is still debated due to the limited theoretical 
approaches at zero temperature. So far rigorous studies of the 


LOFF phase at zero temperature are restricted to its ID struc¬ 
tures including the Fulde-Ferrell (FF) state with a plane-wave 
form A(z) = Ae^'^‘ and the Larkin-Ovhinnikov (LO) state with 
an antipodal-wave form A(z) = 2Acos(2g'z). A recent self- 
consistent treatment of the ID modulation S show that a 
solitonic lattice is formed near the lower critical held, and the 
phase transition to the BCS state is continuous. Near the up¬ 
per critical held the gap function becomes sinusoidal, and the 
transition to the normal state is of first order. 

In addition to these ID structures, there exist a large num¬ 
ber of 3D crystal structures. The general form of a crystal 
structure of the order parameter can be expressed as 

p 

A(r) = ^ Ae^''?*^ ''. (1) 

k=\ 

A specific crystal structure corresponds to a multi-wave con¬ 
figuration determined by the P unit vectors x\.k{k - 1,2,..., P). 
In general, we expect two competing mechanisms: Increas¬ 
ing the number of waves tends to lower the energy, but it may 
also causes higher repulsive interaction energy between differ¬ 
ent wave directions. In a pioneer work, Bowers and Rajagopal 
investigated 23 different crystal structures by using the GL ap¬ 
proach where the grand potential measured with respect 
to the normal state was expanded up to the order O(A^), 

6£l(A) 9 1 /. 1 « 

= PaA^ + -pA^ + -yA® + fKA*") (2) 

with Np being the density of state at the Fermi surface and 
the pair momentum fixed at the optimal value q - \.\9916p. 
Among the structures with y > 0, the favored one seems to be 
the body-centered cubic (BCC) with P - 6 iH. Further, 
it was conjectured that the face-centered cubic (FCC) with 
P - 8 is the preferred structure since its y is negative 
and the largest E3]. For BCC structure, the GL analysis up 
to the order (9(A®) predicts a strong first-order phase transi¬ 
tion at 6 pt - 3.6Ao with the gap parameter A O.SAq in. 
The prediction of a strong first-order phase transition may in¬ 
validate the GL approach itself. On the other hand, by using 
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the quasiclassical equation approach with a Fourier expansion 
for the order parameter, Combescot and Mora [3 pre¬ 
dicted that the BCC-normal transition is of rather weak first 
order: The upper critical field is only about 4% higher 
than 6/^2 with A ^ O.IAq at d/r = d/r*. If this result is reli¬ 
able, it indicates that the higher-order expansions in the GL 
analysis is important for quantitative predictions. To under¬ 
stand this intuitively, let us simply add the eighth-order term 
|A^ to the GL potential ([2]). A detailed analysis of the in¬ 
fluence of a positive 77 on the phase transition is presented in 
Appendix A. We find that with increasing 77 , the first-order 
phase transition becomes weaker and the upper critical field 
d/r, decreases. For 77 —> - 1 - 00 , the phase transition approaches 
second order and 6 fit —> 6 fi 2 . Therefore, to give more precise 
predictions we need to study the higher-order expansions and 
the convergence property of the GL series, or use a different 
way to evaluate the grand potential without assuming a small 
value of A. 

For a specific crystal structure given by O, it is periodic in 
coordinate space. As a result, the eigenvalue equation for the 
fermionic excitation spectrum in this periodic pair potential, 
which is known as the Bogoliubov-de Gennes (BdG) equa¬ 
tion, is in analogy to the Schrodinger equation of quantum par¬ 
ticles in a periodic potential. This indicates that the fermionic 
excitation spectrum has a band structure, which can be solved 
from the BdG equation. The grand potential can be directly 
evaluated once the fermionic excitation spectrum is known 
In this work, we present a solid-state-like calculation 
of the grand potential of the BCC structure. Our numerical 
results show that the phase transition from the BCC state to 
the normal state is of rather weak first order, consistent with 
the work by Combescot and Mora This implies that 

it is quite necessary to evaluate the higher-order terms in the 
GL expansion to improve the quantitative predictions. 


II. THERMODYNAMIC POTENTIAL 


To be specific, we consider a general effective Lagrangian 
for two-flavor quark pairing at high density and at weak cou¬ 
pling. The Lagrangian density is given by 111] 

£ = i//^[id, - e(p) -H fiji/i + £inu (3) 

where i/r = (i/tu, i/id)^ denotes the two-flavor quark field and 
e(p) is the quark dispersion with the momentum operator p = 
-fV. In the momentum representation we have s(p) = jpj. 
The quark chemical potentials are specified by the diagonal 
matrix jl - diag( 7 /u, A<d) in the flavor space, where 


fia - fi + 6fi, fid-fi- 6fi. (4) 

The interaction Lagrangian density which leads to Cooper 
pairing between different flavors can be expressed as 

£mt = g(lAV2lA*)(lA^O-2lA)- (5) 


antiquark degree of freedom because it plays no role at high 
density and at weak coupling. We have also neglected the 
color and spin degrees of freedom, which simply give rise to 
a degenerate factor. 

Color superconductivity is characterized by nonzero ex¬ 
pectation value of the diquark field ^(f, r) = - 2 igi/£cr 2 il/. 
For the purpose of studying inhomogeneous phases, we set 
the expectation value of (p{t,r) to be static but inhomoge¬ 
neous, i.e., r)) = A(r). With the Nambu-Gor’kov spinor 

'y - (if) if)*)^, the mean-field Lagrangian reads 


Xmf = 


id, - e(p) H- p -70-2A(r) \ ^ 

7o-2A*(r) id, + e(p) -ft) 4g 


The order parameters of the BCC and FCC structures can be 
expressed as 

A(r) = 2A [cos (2qx) + cos i2qy) + cos (2qz)\ (7) 


and 


A(r) = 8 Acos 



cos 




( 8 ) 


respectively. Therefore, we consider a 3D periodic structure 
where the unit cell is spanned by three linearly independent 
vectors ai = atx, St 2 - aty, and 33 = with a - nfq 
for BCC and a - '^nfq for FCC. The order parameter is 
periodic in space, i.e., A(r) = A(r + a,). It can be decomposed 
into a discrete set of Fourier components. 


CO 

A(r) = 2AGe'‘^‘'= ^ .(9) 

G l,m,n=—oo 


where the reciprocal lattice vector G reads 

G = G[i„,n] = — (l^x + mCy + ne^ , l,m,n eZ. (10) 
The Fourier component Ag = A[;,„„] can be evaluated as 


- A [ -I- ^7,-1) + difi { 6 ,„d + 




( 11 ) 


and 


Ag - A (did -I- (5;,-i) ((5m ,1 -•- (5m,- 1 ) ((5n,i -I- d„-i) (12) 

for BCC and FCC structures, respectively. 

Then we consider a finite system in a cubic box defined as 
x,y, z 6 [-L/2, L/2] with the length L = Ac 7 . For convenience 
we impose the periodic boundary condition. The thermody¬ 
namic limit is reached by setting N —) oa. Using the momen¬ 
tum representation, we have the Fourier transformation 

'F(t, r) = Yj (13) 

v.p 

Here V is the volume of the system, aiy - (2v+ l) 7 rT(v e Z) is 
the fermion Matsubara frequency, and the quantized momen¬ 
tum p is given by 


where g is the coupling constant and 0-2 is the second Pauli 
matrix in the flavor space. Notice that we have neglected the 


p = 


(14) 
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with l,m,n e Z. The partition function of the system is given Substituting this expansion into the BdG equation, for a given 
by k we obtain a matrix equation 


J [c/'T][c/'T^]e-'^ (15) 

fl/r r -1 

with the Euclidean action S - - Jy “ The grand 

potential per volume reads 

Q^-^lnZ. (16) 

In the mean-field approximation, the action S is quadratic. 
Therefore, the partition function Z and grand potential Q can 
be evaluated. Using the Fourier expansions for T* and A, we 
obtain the mean-field action 

‘5mf = j: 2r ^ (*^vbp,p' - 'Hp,p')'f'vp/,(17) 


2 'KG,G'(k)0G'(k) = £.i(k),^G(k), (23) 

G' 

where the matrix "NG G'ik) is given by 

(|k H-G| -/l)bG,G' /cTaAG-G' \ 

-/cr2A^,_e -(|k + G| - fi)6G,G' ) ■ ^ ^ 

This shows that, for a given k-point in the BZ, we can solve 
the eigenvalue spectrum {^^(k)) by diagonalizing the matrix 
'7/G,G'(k). Without loss of generality, the BZ can be chosen 
as kx, ky, k^ 6 [—71/a, nja]. For a quantized volume V contain¬ 
ing unit cells, we have allowed momenta k in the BZ. 
Accordingly, the grand potential is now given by 

■* keBZ A 


where the effective Hamiltonian matrix TTp p/ reads 

yf, _( (IP|-A)^P,P' 'O'Z Zg ^G^G,P-P' \ /los, 

P-P' \-icr2i:G^h^G,p'^P -(IPl-A)^P,P' /■ 

The effective Hamiltonian TYp p- is a huge matrix in Nambu- 
Gor’kov, flavor, and (discrete) momentum spaces. It is Her- 
mitian and can in principle be diagonalized. Assuming that 
the eigenvalues of is denoted by E^, we can formally 
express the grand potential as 


where the function 'WiE) - | -t- T ln(l -H The summa¬ 

tion over G can be worked out as Zg I^gP = FA^. 

In practice, diagonalization of the matrix TYp p' is infeasi¬ 
ble. However, Ei can be brought into a block-diagonal form 
with independent blocks in the momentum space accord¬ 
ing to the famous Bloch theorem j^ . To understand this, we 
consider the eigenvalue equation for the fermionic excitation 
spectrum in the coordinate space, which is known as the BdG 
equation. For our system, the BdG equation reads 


I e(-/V) - A fcr 2 A(r) 

( -/o-2A*(r) -e(-/V)-i-A 


0 i(r) = £i</'i(r). 


( 20 ) 


According to the Bloch theorem, the solution of the eigen¬ 
function (f>A{r) takes the form of the so-called Bloch function. 
We have 


<p,i(r) = e''‘V 2 k(r), (21) 


In the thermodynamic limit A —> oo, the summation ^ ZkeBZ 
is replaced by an integral over the BZ. 

The Hamiltonian matrix TYG.GAk) can be further simplified 
to lower the matrix size. After a proper rearrangement of the 
eigenvector (pG, we find that EP can be decomposed into two 
blocks. We have EP = EPty^sn ©'^-A.-iSp- The blocks can be 
expressed as EP^Sfi - -Sfi I where I is the identity matrix 
and the matrix EPa is given by 


CHa)g,G' 


(|k -H G| - /i)bG,G' Ag-G' 

^G'-G “(1*^ -H G| - /i)bG,G' 


■ (26) 


The eigenvalues of EPa^si^ do not depend on the sign of A. 
Moreover, replacing the Sjj. by -S/a amounts to a replacement 
of the eigenvalue spectrum {^^(k)) by {-^^(k)). Therefore, 
the two blocks contribute equally to the grand potential and 
we only need to determine the eigenvalues of EPa,sii- The 
Hamiltonian matrix (|2^ represents the general problem of 
two-species pairing with mismatched Fermi surfaces. In the 
weak coupling limit, the pairing is dominated near the Fermi 
surfaces. Therefore, the physical result should be universal 
in terms of the pairing gap Aq at vanishing mismatch and the 
density of state A/p at the Fermi surface. 

In the following we shall focus on the zero-temperature 
case. The grand potential Q is divergent and hence a proper 
regularization scheme is needed. Since we need to deal with 
the Bloch momentum k + G, the usual three-momentum cut¬ 
off scheme is not appropriate for numerical calcula¬ 

tions. Moreover, we are interested in the grand potential 6Q. 
measured with respect to the normal state. Therefore, we em¬ 
ploy a Pauli-Villars-like regularization scheme, in which 6Q. 
is well-defined S. The “renormalized” grand potential is 
given by IH] 


where k is the momentum in the Brillouin zone (BZ) and the 
function <^ 2 k(r) has the same periodicity as the order parame¬ 
ter A(r). We therefore have the similar Fourier expansion 

<^,ik(r) = ^<^G(k)F‘^--. ( 22 ) 

G 


where 


(5G(A, q) - Q(A, q) - G(0, q), 


(27) 


phE 1 r (7% ^ j—. 


,,(k) + ;A2 (28) 
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FIG. 1: (Color online) The lower and upper critical fields (upper 
panel) and the size of the stability window ( 5^2 — <5l^i)/Ao (lower 
panel) for the FF state as a function of Aq at /r = 400 MeV. The thin 
lines denote results in the weak-coupling limit. The blue solid and 
red dashed lines correspond to A = 400 MeV and A = 800 MeV, 
respectively. 


with cq - C 2 — I and ci = -2. Here {^^((k)) denotes the eigen¬ 
value spectrum of The coupling constant g can be fixed 

by the BCS gap Aq at d/i - 0. We expect that at weak coupling 
the physical results depend on the cutoff A only through the 
BCS gap Ao lH^. In Fig. [1] we show the stability window for 
the FF state as a function of Aq. In the weak coupling limit, 
the critical fields depend only on Aq. For accuracy reason 1^ . 
we shall choose Aq ~ 100 MeV at fi - 400 MeV, which cor¬ 
responds to the realistic value of Aq at moderate density 12^ . 
Since the size of the FF window { 6 fj .2 - depends very 

weakly on Aq and A, we can use the upper critical field 6 fi 2 
obtained in Fig. [T]to “calibrate” dyu and appropriately extrap¬ 
olate the results to the weak coupling limit. 


value equation can be rewritten as 

Xi('^A)G.G'(k)<^G'(k) = [£4(k) + 6 IJ] 0G(k), 
G' 


where the Hamiltonian matrix ('KA)G,G'(k) reads 




(29) 


(30) 


Here fp = Ipl - and we have used the fact A^ = A_g- The 
eigenstate <pG includes two components uq and vq as usual in 
the BCS theory. We have 


We notice that 5/i can be absorbed into the eigenvalues. It is 
easy to prove that the eigenvalues of "Ka do not depend on the 
sign of A. Moreover, if s is an eigenvalue of 'T/a, -s must 
be another eigenvalue. Therefore, replacing the dfi by -d// 
amounts to a replacement of the eigenvalue spectrum {/(^((k)) 
by {-£^(k)). 

However, the matrix 'T/a has infinite dimensions because 
the integers l,m,n run from -oo to -Hoo. Therefore, we have 
to make a truncation in order to perform a calculation. It is 
natural to make a symmetrical truncation, i.e., 

-D<l,m,n<D, (D e Z+). (32) 


For sufficiently large D, the contribution from the high-energy 
bands becomes vanishingly small and the grand potential 60. 
converges to its precise value. After making this truncation, 
the matrix equation can be expressed as 


H 


u 

V 


Hn HiaWu 
H21 H22 / \ y 




(33) 


where u and v are {2D + l)^-dimensional vectors and H,y are 
{2D + 1)^ X {2D + 1)^ matrices. The matrix elements of H,y 
can be formally expressed as 


H 

H 


[l,m,n],[l' ,m' ,n'] 
11 

[/,m, 

12 




‘22 

H [l,m,n],[l' ,m' ,n'] 

, 


21 


= A 


[l-l',m-m',n~n ']» 


n,n' ? 

(34) 


where 




2kI 


+ ky H- 


2 mn 


-i-1 k- -I- 


27m\ 
a I 


■M35) 


III. MATRIX STRUCTURE 

For a given k-point in the BZ, we can diagonalize the 
Hamiltonian matrix 'T/a, 5 p to obtain its eigenvalue spectrum 
{/(’^(k)). The choice of the k-points in the BZ should be dense 
enough to achieve the thermodynamic limit . The eigen- 


Here the matrix index [/, m, n] corresponds to the reciprocal 
lattice vector = {2nla){hx + mCy + ne^). It shows that 

the blocks Hu and H 22 are diagonal. The off-diagonal blocks 
Hi 2 and H 21 carry the information of the order parameter A 
and characterize the crystal structure. 

For a specific value of D, we can write down the explicit 
form of the vectors u and u and the matrices H,j. Here we use 
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D = 1 as an example. The vectors u and v are 27-dimensional 
can be expressed as 


where size of the blocks H+ and H_ are both {2D -(-1)^. The 
eigenvector (p is now defined as 


“[-1,-1] ' 


' *^[-1,-1]' 

“[-1,0] 


^<[-1,0] 

“[-1,1] 


f^[-l,l] 

“[0,-1] 


*^[0,-1] 

“[0,0] 

, V - 

*^[0,0] 

“[0,1] 


*^[0,1] 

“[1,-1] 


t^[l,-l] 

“[1,0] 


^[1,0] 

“[1,1] 


V r»[i,i] 


where and V[im\ are defined as 


(36) 



' “[l,m,-l] ' 


' l^[;,m,-l] ' 










In this representation, the off diagonal blocks H 12 and H 21 are 
given by 



Ai 

^2 

0 

A 2 

0 

0 

0 

0 

0 


^2 

Ai A 2 

0 

A 2 

0 

0 

0 

0 


0 

A 2 

^1 

0 

0 

A 2 

0 

0 

0 


^2 

0 

0 

Ai 

A 2 

0 

A 2 

0 

0 

Hi2 = 

0 

A 2 

0 

A 2 

Ai 

A 2 

0 

A 2 

0 


0 

0 

^2 

0 

A 2 

Ai 

0 

0 

A 2 


0 

0 

0 

A 2 

0 

0 

Ai 

A 2 

0 


0 

0 

0 

0 

A 2 

0 

A 2 

Ai 

A 2 


0 

0 

0 

0 

0 

A 2 

0 

A 2 

Ai 

for BCC structure 

and 









0 

0 

0 

0 

Ai 

0 

0 

0 

0 


0 

0 

0 

Ai 

0 

Ai 

0 

0 

0 


0 

0 

0 

0 

Ai 

0 

0 

0 

0 


0 

Ai 

0 

0 

0 

0 

0 

Ai 

0 

Hi2 = 

Ai 

0 

^1 

0 

0 

0 

Ai 

0 

Ai 


0 

Ai 

0 

0 

0 

0 

0 

Ai 

0 


0 

0 

0 

0 

Ai 

0 

0 

0 

0 


0 

0 

0 

Ai 

0 

Ai 

0 

0 

0 


0 

0 

0 

0 

Ai 

0 

0 

0 

0 

for FCC structure 

respectively 

. Here the blocks Aj 

defined as 










/ 

0 

A 0 ) 





{^ 

0 

0 ) 

Ai = 

A 

0 A 



^2 - 

0 

A 

0 

\ 

0 

A 0 , 





lo 

0 

A j 


(38) 


(39) 


(40) 


In principle, the eigenvalue spectrum {^^(k)) can be obtained 
by diagonalizing the matrix H with a size 2{2D H- 1)^. 

We notice that the matrix size 2(2D-i-l)^ grows dramatically 
with increasing cutoff D. Therefore, for realistic diagonaliza- 
tion, it is better to reduce the size of the matrix. Here we find 
that, with a proper rearrangement of the basis (p or after a sim¬ 
ilarity transformation, the matrix H becomes block diagonal. 
We have 


H 


H+ 0 \ 
0 


(41) 


0 = 


U- 


(42) 


For D - I, the 27-dimensional vectors cp^ and (p^ are given by 


' !2[-l _1-1] ' 


' M[-l-1-1] ' 

“[-1,-1,0] 


12[-1,-1,0] 

12[-1,-1,1] 


“[-1,-1,1] 

“[-1,0,-1] 


1^[-1,0,-1] 

12[-1,0,0] 


“[-1,0,0] 

“[-1,0,1] 


1^[-1,0,1] 

12[-1,1,-1] 


“[-1,1,-1] 

“[-1,1,0] 


1^[-1,1,0] 

^^[-1,1,1] 


“[-1,1,1] 

“[0,-l,-l] 


I^IO,-!,-!] 

12[0,-1,0] 


“[0,-1,0] 

“[0,-1,1] 


1^[0,-1,1] 

y[0,0,-l] 


“[0,0,-l] 

“[0,0,0] 

, 0- = 

t^[0,0,0] 

1^[0,0,1] 


“[0,0,1] 

“[0,1,-1] 


1^[0,1,-1] 

1^[0,1,0] 


“[0,1,0] 

“[0,1,1] 


1^[0,1,1] 

'^[1,-1.-1] 


“[1,-1,-1] 

“[1,-1,0] 


1^[1,-1,0] 

1^[1,-1,1] 


“[1,-1,1] 

“[1,0,-1] 


1^[1,0,-1] 

1^[1,0,0] 


“[1,0,0] 

“[1,0,1] 


1^[1,0,1] 

J^[l,l,-1] 


“[1,1,-1] 

“[1,1,0] 


1^[1,1,0] 

, t^[l,l,l] , 


. “[1,1,1] > 


which is just a proper rearrangement of the original basis 
given by (O. The blocks H+ and H_ are given by 


Hi = +Ho + Hi 2 , (44) 

where H 12 is given by ( |34| ) or and for D = 1. Ho is 
a diagonal matrix containing the kinetic energies We 

have 






For D - 1, we obtain 


(45) 


Ho - diag(-^[-i_i_i],^[-i_i,o],-^[-1-1,1], ••• , 

^[0,0,0],-•• ,-^[1,1,-1],^[1,1,0]>-^[1,1,1])- (46) 


It is easy to show that the eigenvalue spectra of H+ and H_ 
are dependent: If the eigenvalue spectrum of H+ is given by 
{£ 2 (k)), the eigenvalue spectrum of H_ reads {-e,)(k)). There¬ 
fore, we only need to diagonalize the matrix H+ or H_ which 
has a size {2D -H 1)^. Once the eigenvalue spectrum of H+ is 
known, the eigenvalue spectrum of Hamiltonian matrix 
is given by 

{£2(k)} = {e2(k) - (5//) U {-e2(k) - 6^i}. (47) 
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FIG. 2: (Color online) The potential curves (511(A) of the BCC structure at the optimal pair momenta for various values of S/j/Aq. The grand 
potential is scaled by a constant Qq = 2.5 x 10*(MeV)'^. The red dots show the data obtained from our numerical calculation. 


Therefore, we can in principle calculate the potential land¬ 
scape AQ(A, q). The solution (A, q) of a specific crystal struc¬ 
ture corresponds to the global minimum of the potential land¬ 
scape. 

IV. COMPUTATION AND RESULTS 

To achieve satisfying convergence we normally need a large 
cutoff D. However, the matrix size (2D -H 1)^ and hence the 
computing time and cost grow dramatically with increasing D. 
The cutoff D needed for convergence can be roughly estimated 
from the maximum momentum in the matrix, 

k^^^(2D+\)- (48) 

a 

The value of A:niax can be estimated from the LO state. The 
calculation of the LO state is much easier than 3D structures 
because the matrix size becomes 2D + 1. The details of the 
calculation of the LO state are presented in Appendix B. For 
Ao ~ 100 MeV we need - 5GeV ll^ . Since we are 
interested in the region byu/Ao e [0.7,0.8] and the optimal 
pair momentum is ^ ~ bp, we estimate D ~ 35 for BCC 
and D ~ 60 for FCC. These huge matrix sizes are beyond the 
capability of our current computing facilities. On the other 
hand, even though a supercomputer may be able to diagonal¬ 
ize these huge matrices, the computing time and cost are still 
enormous, which makes the calculation infeasible. 



(6^-6112)1^0 


FIG. 3: (Color online) Comparison of the grand potentials of various 
phases: BCS (black solid), FF (blue dotted), LO (green dash-dotted), 
and BCC (red dashed). The horizontal axis has been “calibrated” by 
using the quantity (bp - 6 / 12)1 Aq. 


Since we are interested in the grand potential 6D. rather than 
the band structure (the eigenvalues), we can neglect a small 
amount of the off-diagonal couplings A in the matrix H+. By 
doing so, the huge matrix H+ can be decomposed into a num¬ 
ber of blocks with size (2d - 1 -1)^. For symmetry reason, we set 
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the centers of these blocks at the reciprocal lattice vectors 


V. SUMMARY 


- {2d + 1)—(hiCj + riyty + njCj) (49) 

with nx,iiy,n^ e Z. With increasing d, the grand potential 
converges to the result from exact diagonalization. Good con¬ 
vergence is normally reached at some value d - do. The de¬ 
tails of our computational scheme are presented in Appendix 
C. If the block size (2do -H 1)^ is within our computing capa¬ 
bility, the calculation becomes feasible. Fortunately, we find 
that this computational scheme works for the BCC structure. 
At present, we are not able to perform a calculation for the 
FCC structure, since the value of do needed for convergence 
is much larger. Note that the computing cost is still very large 
even though we have employed this effective computational 
scheme. 

We have performed calculations of the BCC structure for 
A„ = 60,80,100 MeV iH at /r = 400 MeV il. For differ¬ 
ent values of Aq, the results are almost the same in terms of 
the quantity (dfj. - dfi 2 )l^o- Therefore, we anticipate that our 
results can be appropriately extrapolated to the weak coupling 
limit Ao ^ 0. In the following, we shall present the result 
for Ao = 100 MeV. For a given value of b/z/Ao, we calculate 
the potential curve bf2(A) at various values of q and search for 
the optimal pair momentum and the minimum of the potential 
landscape. The potential curves at the optimal pair momenta 
for several values of b/z/Ao are shown in Fig. |2 With increas¬ 
ing value of b/z/Ao, the potential minimum gets shallower. At 
a critical value dfit — 6 fi 2 - 0.03Ao, the potential minimum 
approaches zero and a first-order phase transition to the nor¬ 
mal state occurs. The comparison of the grand potentials of 
various phases are shown in Fig. |3] For the LO state, its phase 
transition to the normal state occur almost at the same point 
as the FF state, b;U 2 - O.SAq. At b/z = b/Z 2 , the grand poten¬ 
tial of the BCC structure is negative, which indicates that the 
BCC structure is energetically favored around the FF-normal 
transition point. Well below the FF-normal transition point, 
the BCC state has higher grand potential than the LO state 
and hence is not favored. Near the BCS-LO transition, the 
solitionic state becomes favored S. However, this does not 
change our qualitative conclusion. 

Our result is qualitatively consistent with the GL analy¬ 
sis d. However, the quantitative difference is significant: 
The GL analysis predicts a strong first-order phase transition 
and a large upper critical field lUOl] . while our result shows 
a weak first-order phase transition at which A ^ 0.1 Aq. On 
the other hand, our result is quantitatively compatible with the 
quasiclassical equation approach where it shows that 

the BCC structure is preferred in a narrow window around 
byU = 6 fi 2 at zero temperature d. Therefore, the GL analy¬ 
sis up to the order (9(A®) may not be quantitatively sufficient. 
We notice that the LO state already shows the limitation of the 
GL analysis: While the GL analysis predicts a second-order 
phase transition, exact calculation shows a first-order phase 
transition (see also Appendix B). In the future, it is neces¬ 
sary to study the higher-order expansions and the convergence 
property of the GL series, which would help to quantitatively 
improve the GL predictions. 


In summary, we have performed an solid-state-like calcu¬ 
lation of the ground-state energy of a 3D structure in crys¬ 
talline color superconductivity. We proposed a computational 
scheme to overcome the difficulties in diagonalizing matrices 
of huge sizes. Our numerical results show that the BCC struc¬ 
ture is preferred in a small window around the conventional 
FF-normal phase transition point, which indicates that the 
higher-order terms in the GL approach are rather important. In 
the future it would be possible to perform a calculation for the 
FCC structure with stronger computing facilities and/or with 
better method of matrix diagonalization. This solid-state-like 
approach can also be applied to study the crystalline structures 
of the three-flavor color-superconducting quark matter El 
and the inhomogeneous chiral condensate ESl- 
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Appendix A: Ginzburg-Landau Theory: Importance of Higher 
Order Expansions 

In the Ginzburg-Landau (GL) theory of crystalline color 
superconductors at zero temperature and at weak coupling, 
the grand potential measured with respect to the normal state, 
bQ = Q - Gn, is expanded as El 

^ = PaA^ V + + ^(A'O), (A1) 

A/^f 2 3 4 

where TVp is the density of state at the Fermi surface. The 
coefficient a is universal for all crystal structures and is given 

by El 


b/z H- b|z 1 
a - -1 -In-In ■ 


A2 


2 q ^ - b/z 2 4{q^ - 


(A2) 


Let us consider the vicinity of the conventional LOFF-normal 
transition point b/z = b/Z 2 , where we have 

— = 0.7544, — = 1.1997. (A3) 

Ao b/Z2 

At the pair momentum q - 1.1997b/z, we obtain 
b/z , b/Z 2 , b/z 

a = In -In —— = In-. (A4) 

Ao Ao b/Z2 
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FIG. 4: The GL potential curves of the BCC structure for different values of rj alSfi = 6^2 = 0.754Ao. 




A/Ao A/Ao 

FIG. 5: The GL potential curves of the BCC structure for different values of j) at the first-order phase transition point 6 ^ = 6 ^,. 


For convenience, we make the GL potential dimensionless by 
using the variables 5Q = 6 D.I{NQd^^), A - A/dfi, ]3 = 
y - and 77 = We have 

dO = PaA^ + H- iyA'’ H- ^fjA^ + 0(A^^). (A5) 

The GL coefficients p and y for a number of crystalline 
structures were first calculated by Bowers and Rajagopal ifl^ . 
The predictions for the nature of the phase transitions were 
normally based on the GL potential up to the sixth order (A®). 
To our knowledge, the higher order GL coefficients have never 
been calculated. Here we show that the higher-order GL ex¬ 
pansions are important for the prediction of the phase transi¬ 
tion. To be specific, let us consider the BCC structure. Its GL 
coefficients p and y have been evaluated as d 

^6 =-31.466, 7=19.711. (A 6 ) 

Since p <0, the phase transition to the normal state should be 
of first order. If we employ the GL potential up to the sixth 
order, we predict a strong first-order phase transition at 6 fi - 
6 lit — 3.625Ao. Let us turn on the eighth-order term and study 
how the size of the coefficient ij influences the quantitative 
prediction of the phase transition. In Fig. |4] we show the 


GL potential curves for two different values of 77 at = 6112 - 
For vanishing 77, the potential curve develops a deep minimum 
^Qmin - -13.4 at A = O. 95 A 0 , which indicates a strong first- 
order phase transition at dfi = dyu, » 6112 - However, for a 
large value 77 = 1000, we find a shallow minimum dGmin - 
-0.21 located at A = 0.31Ao. In Fig. |5] we show the GL 
potential curves at the first-order phase transition point dn - 
5Ilf For ?7 = 0 we find a strong first-order phase transition 
at 611 - 6 n, = 3.625Ao, where the minima located at A = 0 
and A = 0.83Ao become degenerate. For fj - 1000, however, 
we observe a much weaker first-order phase transition at bp = 
bp, = 0.95 lAo, where the degenerate minima are located at 
A = 0 and A = 0.28Ao. These results clearly show that, for 
larger 77, the first-order phase transition becomes weaker. For 
77 ^ -1-00, we expect that bp, ^ bp 2 = 0.754Ao. On the other 
hand, if fj is small or even negative, then the next order A*° 
would become important. 


Appendix B: Calculation of the LO State 

The order parameter of the LO state is given by 
A(z) = 2Acos(2g'z). 


(Bl) 
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It is periodic along the z direction with the periodicity a - 
711q. So it can be decomposed into a discrete set of Fourier 
components, 


OO 

A(z) = ^ (B2) 

n=-oo 


The Fourier component A„ is given by 

A„ = - r *A(z)e-2«'«^-A(d„.i+5„,_i). (B3) 

a Jo 


The matrices H,y 

can be explicitly written 

as 





f-2 

0 

0 

0 





0 

4- 

0 

0 

Hn 

= -H 22 - 



0 

0 

^0 

0 





0 

0 

0 





\ 

0 

0 

0 

0 



/ 

0 

A 

0 

0 0 

\ 




A 

0 

A 

0 0 


Hi2 

= H21 = 


0 

A 

0 

A 0 





0 

0 

A 

0 A 




\ 

0 

0 

0 

A 0 



0 f 
0 
0 
0 

. 


(BIO) 


The matrix equation takes a similar form as the 3D structures. 
We have 

^CHA)„.„Tk)0„Tk) = [Ej(k) + 6jj] <p„(k), (B4) 

n' 


where the Hamiltonian matrix ('KA)H,H'(k) reads 


The eigenvalue spectrum {£’,i(k)) can be obtained by diago¬ 
nalizing the matrix H with a size 2{2D -i- 1). With a proper 
rearrangement of the basis or a similarity transformation, 
we have 


H 


H+ 0 \ 
0 H_ j’ 


(Bll) 




^n^n,n’ ^n—n' 

A/2— ^n^n,n' 


with 


(B5) 


where the sizes of H+ and H_ are both 2D + 1. The basis <p is 
now defined as 


0 = 


/ 

U- 


(B12) 


= ^jkl+(k, + 2nqf-ii. (B6) 

We notice that only the motion in the z direction becomes 
quantized. The BZ for can be defined as -tt/o < < Jrla 

or —q<k^< q. The eigenstate includes two components 
M„ and v„. We have 


For D -2, the 5-dimensional vectors 0+ and are given by 


' M-2 ' 


' V-2 ' 

V-l 


M_1 

Uq 


Vo 

Vi 


Ml 

U2 y 


. V2 . 


0«(k) = 


M«(k) \ 

! 2 „(k )) ■ 


(B7) 


The blocks H+ and H_ can be expressed as 


Hi = +H„ + Hi2 . 


(B14) 


If e is an eigenvalue of TTa, -e must be another eigenvalue. 
Therefore, replacing the 6/j by -b/r amounts to a replacement 
of the eigenvalue spectrum {^^(k)) by {-^^(k)). 

After a truncation -D < n < D, we obtain a finite matrix 
equation 


H 


u 

V 


H„ 

H21 H22 / y *2 




(B8) 


where u and v are {2D + l)-dimensional vectors and H,j are 
{2D -H 1) X {2D -Hi) matrices. For a specific value of D, we can 
write down the explicit form of the vectors u and v and the 
matrices H,j. Here we use D = 2 as an example. The vectors 
u and V are 5-dimensional can be expressed as 


' U-2 ' 


' V-2 ' 

U-\ 


V-i 

Uo 

, V = 

VQ 

Ml 


Vl 

1 M2 . 


. V2 . 


Ho is a diagonal matrix containing the kinetic energies. We 
have 


(Ho),,„- = {-mn6„,„. (B15) 


The eigenvalue spectra of H+ and H_ are dependent: If the 
eigenvalue spectrum of H+ is given by {e,i(k)), the eigenvalue 
spectrum of H_ reads {-e,i(k)). Therefore, we need only 
to diagonalize the matrix H+ or H_ which has a dimension 
2D + 1. Once the eigenvalue spectrum of H+ is known, the 
eigenvalue spectrum of Hamiltonian matrix is given by 
{^(k)) = {e2(k) - bju) U {-ei(k) - 

The thermodynamic potential of the LO state at zero tem¬ 
perature can be expressed as 




Similar Pauli-Villas-like regularization scheme should be ap¬ 
plied finally. In Fig. |6] (a), we show the grand potential of 
the LO state for Aq = 80 MeV. The grand potential for the 
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FIG. 6: (Color Online) (a) Comparison of the grand potentials of the LO state and the self-consistent ID modulation for Aq = 80 MeV. (b) The 
potential curve of the LO state at d/r = 0.775Ao and at the optimal pair momentum q = 0.9Ao or q = 1.1613d/j. 


self-consistent ID modulation for Aq = 80 MeV was also re¬ 
ported in ii. We find that the results for the LO state and the 
self-consistent ID modulation agrees with each other near the 
phase transition to the normal state. Near the BCS-LO transi¬ 
tion point, the self-consistent ID modulation has lower grand 
potential than the LO state. It was shown in ll^ that the self- 
consistent ID modulation forms a soliton lattice structure near 
the lower critical field, which lowers the grand potential of the 
system. Near the upper critical field the gap function becomes 
sinusoidal, and therefore the grand potentials of the LO state 
and the ID modulation agree with each other. We notice that 
the phase transition from the LO state to the normal state is 
of first order, which is in contradiction to the prediction from 
the GL analysis. To understand the reason, we show in Fig. |6] 
(b) the potential curve at A// = 0.775Ao and at the optimal pair 
momentum^' = L1613d/i. We find that the potential curve has 
two minima: a shallow minimum at A ^ 0.12Ao and a deep 
minimum at A ^ 0.44Ao. Obviously, the shallow minimum is 
responsible for the GL theory which predicts a second-order 
phase transition. However, the deep global minimum, which 
cannot be captured by the GL theory up to the order A®, is 
responsible for the real first-order phase transition. Therefore, 
the LO state already shows the importance of the higher-order 
expansions in the GL theory. 


Appendix C: Calculation of the Grand Potential: Small Block 
Method 

The key problem in the numerical calculation is the diago- 
nalization of the matrix H+ or H_ and obtaining all the eigen¬ 
values. For BCC and FCC structures, we use a symmetrical 
truncation -D < I, m, n < D with a large cutoff D e Z^. How¬ 
ever, the matrix size grows dramatically with increasing cutoff 
D, which makes the calculation infeasible because of not only 
the computing capability of current computing facilities but 
also the computing time and cost. Notice that we need to di¬ 
agonalize the matrix H+ for various values of the momentum 
k in the BZ, the gap parameter A, and the pair momentum q. 

We first estimate the size of D needed for the convergence 


of the grand potential AQ. The matrix size {2D H-1 and hence 
the computing time and cost grow dramatically with increas¬ 
ing D. The cutoff D is related to the maximum momentum 
kmax in each direction {x, y, and z). We have 

k„ax = (2D+l)-. (Cl) 

a 

This maximum momentum can be roughly estimated from the 
calculation of the LO state. For the LO state, the matrix size 
becomes 2D -i- 1 and exact diagonalization is possible. The 
regime of dfi we are interested in is dfilAo e [0.7 - 0.8] and 
the optimal pair momentum is located at q ^ dfi. From the 
calculation of the LO state at Aq ~ 100 MeV, we find that 
kmax must reach at least 5GeV for convergence. Notice that 
we have k^^Rx - {2-D + \)q for BCC and VSkmax = (2D + 1)^' 
for FCC. Therefore, the cutoff D for BCC can be estimated 
as D ~ 35, which corresponds to a matrix size ~ 3 x 10^. 
For FCC, the cutoff is even larger because of the factor V3. 
We have D ~ 60 for FCC, which corresponds to a matrix size 
~ 1.5 X 10®. Notice that this is only a naive estimation. In 
practice, the cutoff needed for convergence may be smaller or 
larger. Exact diagonalization of such huge matrices to obtain 
all the eigenvalues are impossible with our current computing 
facility. 

We therefore need a feasible scheme to evaluate the grand 
potential AQ. Notice that decreasing the value of Aq does not 
reduce the size of the matrices. In this case, even though k^^x 
becomes smaller, the pair momentum q also gets smaller. Let 
us call an off-diagonal element A in H+ or H_ a “coupling”. 
Because our goal is to evaluate the grand potential AQ rather 
than to know exactly all the band dispersions (eigenvalues), 
we may neglect a small amount of couplings to lower the size 
of the matrices. By neglecting this small amount of couplings, 
the huge matrix H+ becomes block diagonal with each block 
having a much smaller size. In general, we expect that the 
omission of a small amount of couplings A induces only a per¬ 
turbation to the grand potential AQ. We shall call this scheme 
small block method (SBM). 

To be specific, the size of the small blocks in our calculation 
is {2d -1-1)^ with d e 1A. In general, we have d < D. For 
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FIG. 7: Comparison of the grand potentials calculated from the exact diagonalization and from the small block method, (a) The relative error 
R for the LO state at S/j-IAq = 0.77 and = 1.16 with D = 50 and d = 20. (b) The grand potential for the BCC state at 6^1 Ao = 0.75, 
q/6fi = 1.07, and A/Ao = 0.167 as a function of d. 


symmetry reason, we require that the centers of these blocks 
are located at the reciprocal lattice vectors 

In 

= (2t7 + + riy^y + n^^-) (C2) 

with tlx, ny, til s This scheme makes the SBM feasible even 
though (2D + 1)^ is not divisible by (2d +1)^. In practice, 
we first choose a large cutoff D which is sufficient for conver¬ 
gence. By increasing the value of d, we find that the grand 
potential 60. finally converges. In practice, if the grand poten¬ 
tials evaluated at several values of d, i.e., do - k, do - k + 1, ... 
, and do (k 6 Z^), are very close to each other, we identify that 


the grand potential converges to its precise value from exact 
diagonalization. At the converging value d = do, the block 
size (2do + 1)^ is normally much smaller than the total size 
(2D -H 1)^. This scheme makes the calculation feasible and 
also saves a lot of computing time and cost. 

The matrices for the 3D structures are huge and cannot be 
written down here. For the sake of simplicity, let us use the 
LO state as a toy example for the SBM. In this case, the matrix 
size and the block size are 2D + 1 and 2d + 1, respectively. 
The centers of the blocks are located at the reciprocal lattice 
vectors (2d 4 - l)2q'n,ej with e Z. For D = 10, the matrix 
H+ reads 


' e+io A 0 0 0 

A s+9 A 0 0 

0 A e+8 A 0 

0 0 A e+7 A 

0 0 0 A s+6 

0 0 0 0 A 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

, 0 0 0 0 0 


0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

A 0 0 0 0 0 

s +5 A 0 0 0 0 

A e+4 A 0 0 0 

0 A s +3 A 0 0 

0 0 A e+2 A 0 

0 0 0 A A 

0 0 0 0 A So 

0 0 0 0 0 A 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 

0 0 0 0 0 0 


0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

A 0 0 0 0 

e_i A 0 0 0 

A e_2 A 0 0 

0 A £_3 A 0 

0 0 A e_4 A 

0 0 0 A e_5 

0 0 0 0 A 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 


0 0 0 0 Of 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

A 0 0 0 0 

e_6 A 0 0 0 

A e_7 A 0 0 

0 A A 0 

0 0 A e_9 A 

0 0 0 A e_io , 


(C3) 
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where e„ 

= (-1)" 


+ {k~ + 2nq)^ - 


. If we take d = 

2, we neglect the couplings A in red 

approximated as 
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In this case, the matrix H+ is 


0 0 0 1 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

£-8 A 0 

A e_g A 

0 A e_io, 


(C4) 


Therefore, by neglecting a small amount of couplings, we 
have made the large matrix H+ block-diagonal. Notice that 
this is only a toy example for the SBM. In practice, D = 10 
and c/ = 2 is obviously not enough for convergence. 

For the LO state, exact diagonalization of the matrices at 
q ^ 6iu is quite easy because the size of the matrices is 2D + 1. 
We can therefore check the error induced by the SBM. The 
relative error induced by the SBM can be defined as 

R = I^^SBM-dQExI 

where (5f2sBM and (5f2Ex are the grand potentials obtained from 
the SBM and exact diagonalization, respectively. In Fig. |2](a), 
we show a numerical example of the relative error for the LO 
state at i5ju/Ao = 0.77 and q/dj^ = 1.16. In the calculations, 
we use Z) = 50 and d = 20. We find that the relative error 
is very small, generally of order (9(10^^). The slightly larger 
error around A/Ao = 0.5 is due to the fact that dQ itself is 


very small there. For the BCC structure, we are not able to 
check the relative error at q ^ 6fi because it is impossible 
to exactly diagonalize the matrices with a huge size {2D + 
1)^. However, we can check the d dependence of the grand 
potential. For pair momentum around q ^ dji, we choose a 
sufficiently large cutoff D and increase the value of d. We 
evaluate the grand potentials for various values of d (i.e., do - 
k, do - k + 1,, and do). If they are very close to each other, 
we identify that the grand potential converges. Then the grand 
potential 6Q. can be evaluated at d - do. In Fig. |7] we show 
the d dependence of the grand potential of the BCC structure 
at dyu/Ao = 0.75, = 1-07, and A/Aq = 0.167. In the 

calculation we choose D - 50 which is sufficiently large to 
guarantee the convergence at large G. We find that for BCC 
structure, do is normally within the range 10 < c/q < 15, which 
is feasible for a calculation. For FCC structure, we do not find 
a satisfying convergence at these small values of d. 
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